| Plant | Flowers | Date | lon | lat | ele | Year | Month |
|---|---|---|---|---|---|---|---|
| Glossoloma oblongicalyx | 4 | 2015-10-19 | -78.59093 | 0.130838 | 2270 | 2015 | 10 |
| Gasteranthus quitensis | 2 | 2016-10-17 | -78.59770 | 0.120070 | 1940 | 2016 | 10 |
| Kohleria affinis | 1 | 2016-12-13 | -78.59534 | 0.126746 | 2110 | 2016 | 12 |
| Columnea ciliata | 3 | 2014-02-27 | -78.59934 | 0.116682 | 1960 | 2014 | 2 |
| Columnea medicinalis | 1 | 2014-04-23 | -78.59372 | 0.128700 | 2130 | 2014 | 4 |
| Drymonia teuscheri | 3 | 2016-07-28 | -78.59245 | 0.129393 | 2200 | 2016 | 7 |
Let’s start with a traditional randomization test to get everyone comfortable with the essential result. For each random draw, calculate the mean niche overlap in hummingbird usage. Species can only bloom at site which they occur. Mantains plant abundance.
Species elevation models without conditional probability of flowering
## sink("model/occurrence.jags")
## cat("
## model {
##
## for (x in 1:Nobs){
##
## #Observation of a flowering plant
## Y[x] ~ dbern(p[x])
## logit(p[x]) <- alpha[Plant[x]] + beta[Plant[x]] * elevation[x]
##
## #Residuals
## discrepancy[x] <- abs(Y[x] - p[x])
##
## #Assess Model Fit
## Ynew[x] ~ dbern(p[x])
## discrepancy.new[x]<-abs(Ynew[x] - p[x])
## }
##
## #Sum discrepancy
## fit <- sum(discrepancy)/Nobs
## fitnew <- sum(discrepancy.new)/Nobs
##
## #Prediction
##
## for(x in 1:Npreds){
## #predict value
##
## #Observation - probability of flowering
## prediction[x] ~ dbern(p_new[x])
## logit(p_new[x]) <- alpha[NewPlant[x]] + beta[NewPlant[x]] * new_elevation[x]
##
## #predictive error
## pred_error[x] <- abs(Ypred[x] - p_new[x])
## }
##
## #Predictive Error
## fitpred<-sum(pred_error)/Npreds
##
## #Priors
##
## #Species level priors
##
## for (j in 1:Plants){
## #Intercept flowering probability
## alpha[j] ~ dnorm(0,0.386)
## beta[j] ~ dnorm(0,0.386)
## }
##
## }
## ",fill=TRUE)
##
## sink()
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 3600
## Unobserved stochastic nodes: 4416
## Total graph size: 46411
##
## Initializing model
Confirm everything matches up
## # A tibble: 1 x 8
## Plant Site Month Year n ele Index jPlant
## <fct> <dbl> <dbl> <chr> <dbl> <dbl> <chr> <dbl>
## 1 Columnea ciliata 3 3 2018 1 1.78 50 2
## # A tibble: 4 x 7
## Plant Site Month Year n ele Index
## <fct> <dbl> <dbl> <chr> <dbl> <dbl> <int>
## 1 Columnea ciliata 3 3 2014 1 1.77 18
## 2 Columnea ciliata 3 3 2015 1 1.76 66
## 3 Columnea ciliata 3 3 2018 1 1.78 201
## 4 Columnea ciliata 3 3 2019 1 1.79 247
## Draw chain par value parameter jPlant
## 1 1 1 discrepancy[100] 0.6163511 discrepancy 4
## 2 2 1 discrepancy[100] 0.7356973 discrepancy 4
## 3 3 1 discrepancy[100] 0.6173596 discrepancy 4
## 4 4 1 discrepancy[100] 0.6357632 discrepancy 4
## 5 5 1 discrepancy[100] 0.5818248 discrepancy 4
## 6 6 1 discrepancy[100] 0.6252359 discrepancy 4
## Plant Index Site Month Year n ele
## 1 Columnea mastersonii 100 1 12 2016 1 1.364
## 2 Columnea mastersonii 100 1 12 2016 1 1.364
## 3 Columnea mastersonii 100 1 12 2016 1 1.364
## 4 Columnea mastersonii 100 1 12 2016 1 1.364
## 5 Columnea mastersonii 100 1 12 2016 1 1.364
## 6 Columnea mastersonii 100 1 12 2016 1 1.364
Probability of observing a flowering record
Discrepancy
## mean upper lower
## 1 0.6507852 0.7559506 0.5431877
Fit by presence/abs
Probability of flowering conditional on elevation occurrence
## sink("model/occurrence_month.jags")
## cat("
## model {
##
## for (x in 1:Nobs){
##
## #Observation of a flowering plant
## Y[x] ~ dbern(p[x])
## logit(p[x]) <- alpha[Plant[x],Month[x]] + beta[Plant[x]] * elevation[x]
##
## #Residuals
## discrepancy[x] <- abs(Y[x] - p[x])
##
## #Assess Model Fit
## Ynew[x] ~ dbern(p[x])
## discrepancy.new[x]<-abs(Ynew[x] - p[x])
## }
##
## #Sum discrepancy
## fit <- sum(discrepancy)/Nobs
## fitnew <- sum(discrepancy.new)/Nobs
##
## #Prediction
##
## for(x in 1:Npreds){
## #predict value
##
## #Observation - probability of flowering
## prediction[x] ~ dbern(p_new[x])
## logit(p_new[x]) <- alpha[NewPlant[x],NewMonth[x]] + beta[NewPlant[x]] * new_elevation[x]
##
## #predictive error
## pred_error[x] <- abs(Ypred[x] - p_new[x])
## }
##
## #Predictive Error
## fitpred<-sum(pred_error)/Npreds
##
## #Priors
## #Species level priors
## for (j in 1:Plants){
## #effect of elevation
## beta[j] ~ dnorm(0,0.386)
## for(k in 1:Months){
## #Intercept flowering probability by month
## alpha[j,k] ~ dnorm(0,0.386)
## }
##
## }
##
## }
## ",fill=TRUE)
##
## sink()
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 3600
## Unobserved stochastic nodes: 4592
## Total graph size: 51196
##
## Initializing model
Elevation response
month intercepts
Fit by presence/absence
Confirm everything matches up
## # A tibble: 1 x 8
## Plant Site Month Year n ele Index jPlant
## <fct> <dbl> <dbl> <chr> <dbl> <dbl> <chr> <dbl>
## 1 Besleria solanoides 5 9 2014 1 2.14 20 1
## # A tibble: 5 x 7
## Plant Site Month Year n ele Index
## <fct> <dbl> <dbl> <chr> <dbl> <dbl> <int>
## 1 Besleria solanoides 5 9 2014 1 2.14 41
## 2 Besleria solanoides 5 9 2015 1 2.14 88
## 3 Besleria solanoides 5 9 2016 1 2.14 126
## 4 Besleria solanoides 5 9 2017 0 2.23 173
## 5 Besleria solanoides 5 9 2018 0 2.23 225
Probability of observing a flowering record
Discrepancy
## mean upper lower
## 1 0.684529 0.8883212 0.4516712
## sink("model/occurrence_attraction.jags")
## cat("
## model {
##
## for (x in 1:Nobs){
## #Observation of a flowering plant
## Y[x] ~ dbern(p[x])
## logit(p[x]) <- e[Plant[x],Month[x]] + beta[Plant[x]] * elevation[x]
##
## #Residuals
## discrepancy[x] <- abs(Y[x] - p[x])
##
## #Assess Model Fit
## Ynew[x] ~ dbern(p[x])
## discrepancy.new[x]<-abs(Ynew[x] - p[x])
## }
##
## #Sum discrepancy
## fit <- sum(discrepancy)/Nobs
## fitnew <- sum(discrepancy.new)/Nobs
##
## #Prediction
##
## for(x in 1:Npreds){
## #predict value
##
## #Observation - probability of flowering
## prediction[x] ~ dbern(p_new[x])
## logit(p_new[x]) <- e[NewPlant[x],NewMonth[x]] + beta[NewPlant[x]] * new_elevation[x]
##
## #predictive error
## pred_error[x] <- abs(Ypred[x] - p_new[x])
## }
##
## #Predictive Error
## fitpred<-sum(pred_error)/Npreds
##
## #########################
## #autocorrelation in error
## #########################
##
## #For each of month, error covariance among flowering plants
## for(k in 1:Months){
## e[1:Plants,k] ~ dmnorm(zeros,tauC[,])
## }
##
## ##covariance among similar species
## for(i in 1:Plants){
## for(j in 1:Plants){
## C[i,j] = exp(-lambda_cov* D[i,j])
## }
## }
##
## ## Covert variance to precision for each parameter, allow omega to shrink to identity matrix
## vCov = omega*C[,] + (1-omega) * I
## tauC = inverse(vCov*gamma^2)
##
## #Autocorrelation priors
## gamma ~ dunif(0,20)
##
## #Strength of covariance decay
## lambda_cov ~ dunif(0.1,2)
## omega ~ dbeta(1,1)
##
## #Priors
## #Species level priors
## for (j in 1:Plants){
## #effect of elevation
## beta[j] ~ dnorm(0,0.386)
## alpha[j] ~ dnorm(0,0.386)
## }
##
## }
## ",fill=TRUE)
##
## sink()
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 3600
## Unobserved stochastic nodes: 4431
## Total graph size: 52012
##
## Initializing model
Elevation response
month intercepts
## sink("model/occurrence_repulsion.jags")
## cat("
## model {
##
## for (x in 1:Nobs){
## #Observation of a flowering plant
## Y[x] ~ dbern(p[x])
## logit(p[x]) <- e[Plant[x],Month[x]] + beta[Plant[x]] * elevation[x]
##
## #Residuals
## discrepancy[x] <- abs(Y[x] - p[x])
##
## #Assess Model Fit
## Ynew[x] ~ dbern(p[x])
## discrepancy.new[x]<-abs(Ynew[x] - p[x])
## }
##
## #Sum discrepancy
## fit <- sum(discrepancy)/Nobs
## fitnew <- sum(discrepancy.new)/Nobs
##
## #Prediction
##
## for(x in 1:Npreds){
## #predict value
##
## #Observation - probability of flowering
## prediction[x] ~ dbern(p_new[x])
## logit(p_new[x]) <- e[NewPlant[x],NewMonth[x]] + beta[NewPlant[x]] * new_elevation[x]
##
## #predictive error
## pred_error[x] <- abs(Ypred[x] - p_new[x])
## }
##
## #Predictive Error
## fitpred<-sum(pred_error)/Npreds
##
## #########################
## #autocorrelation in error
## #########################
##
## #Error covariance among flowering plants
## for(k in 1:Months){
## e[1:Plants,k] ~ dmnorm(zeros,tauC[,])
## }
##
## ##covariance among similar species
## for(i in 1:Plants){
## for(j in 1:Plants){
## C[i,j] = exp(-lambda_cov*D[i,j])
## }
## }
##
## ## Covert variance to precision for each parameter, allow omega to shrink to identity matrix
## vCov = omega*C[,] + (1-omega) * I
## tauC = vCov*gamma^2
##
## #Autocorrelation priors
## gamma ~ dunif(0,20)
##
## #Strength of covariance decay
## lambda_cov ~ dunif(0.1,2)
## omega ~ dbeta(1,1)
##
## #Priors
## #Species level priors
## for (j in 1:Plants){
## #effect of elevation
## beta[j] ~ dnorm(0,0.386)
## alpha[j] ~ dnorm(0,0.386)
## }
##
## }
## ",fill=TRUE)
##
## sink()
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 3600
## Unobserved stochastic nodes: 4431
## Total graph size: 52011
##
## Initializing model
Glossoloma purpureum
## # A tibble: 4 x 2
## Model p
## <chr> <dbl>
## 1 baseline 0.854
## 2 interaction_attraction 0.599
## 3 interaction_repulsion 0.485
## 4 occurrence_month 0.42
Without baseline
## # A tibble: 3 x 2
## Model p
## <chr> <dbl>
## 1 interaction_attraction 0.599
## 2 interaction_repulsion 0.485
## 3 occurrence_month 0.42
| Model | mean | lower | upper |
|---|---|---|---|
| interaction_repulsion | 0.1894753 | 0.1839611 | 0.1954302 |
| occurrence_month | 0.1846309 | 0.1789345 | 0.1904198 |
| interaction_attraction | 0.1685083 | 0.1628142 | 0.1742645 |
| Model | mean | lower | upper |
|---|---|---|---|
| baseline | 0.1847199 | 0.1783802 | 0.1911261 |
| interaction_repulsion | 0.1769840 | 0.1705929 | 0.1838919 |
| occurrence_month | 0.1734432 | 0.1665216 | 0.1802577 |
| interaction_attraction | 0.1660734 | 0.1592176 | 0.1728697 |
Columnea medicinialis v Columnea strigosa which have strong overlap in visitors
Dint["Columnea medicinalis","Columnea strigosa"]
## [1] 0.3083944
Logit e: Autocorrelation among species
InvLogit e: Autocorrelation among species
## # A tibble: 2 x 4
## # Groups: Model [2]
## Model Var1 Var2 Correlation_E
## <chr> <fct> <fct> <dbl>
## 1 interaction_attraction Columnea medicinal… Columnea strigo… NA
## 2 interaction_repulsion Columnea medicinal… Columnea strigo… NA
## # A tibble: 4 x 4
## # Groups: Model [4]
## Model Var1 Var2 Correlation_P
## <chr> <fct> <fct> <dbl>
## 1 baseline Columnea medicinal… Columnea strigo… 0.991
## 2 interaction_attraction Columnea medicinal… Columnea strigo… 0.963
## 3 interaction_repulsion Columnea medicinal… Columnea strigo… 0.963
## 4 occurrence_month Columnea medicinal… Columnea strigo… 0.774
## # A tibble: 4 x 4
## # Groups: Model [4]
## Model Var1 Var2 Correlation_Ynew
## <chr> <fct> <fct> <dbl>
## 1 baseline Columnea medicina… Columnea strig… NA
## 2 interaction_attracti… Columnea medicina… Columnea strig… 0.768
## 3 interaction_repulsion Columnea medicina… Columnea strig… -0.136
## 4 occurrence_month Columnea medicina… Columnea strig… -0.630
Flowering patterns
Model Fit and Prediction Fit
Comparison of the attraction and repulsion
Overlay observed
Prediction
Predictions from the best model
Supplamental figures